% This file generates the IRF presented in the fifth quadrant (2,2) in 
% Figure 9 in the paper. It needs to be ran from 
% 'FIG_9_Robustness_Output.m'.

global P2 lambda P H
%% Options for the specific results to generate.
graph           =   0;          % = 1 to generate a plot.

csv_file        =   1;          % = 1 if you want to import from CSV.
                                % = 0 if you want to import from XLSX.

use_data        =   1;          % = 1   if baseline data
                                % = 3.1 if output is Industrial Production Index 
                                % = 3.3 if output is hours per capita nonfarm
                                % = 3.5 if output is unemployment rate.  
%% Importing and Preparing the data.

if csv_file == 1
    data            =   csvread('Data_File_Baseline.csv');
    % UM Expectations Median:     
    pilev_data      =   csvread('Data_File_Inflation.csv');
    pilev           =   pilev_data(:,6);     
else
    data            =   xlsread('Data_File.xlsx','Baseline');
    % UM Expectations Median:
    pilev           =   xlsread('Data_File.xlsx','Inflation','F2:F281');
end
% Import the specific variables. 
dz              =   data(:,3);
dy              =   data(:,4);
zlb             =   data(:,5);
rec             =   data(:,6);
zlev            =   data(:,7);
time            =   data(:,9);
time22          =   time.^2;
[T,n]           =   size(dy);

if use_data == 1             % GDP
    ylev        =   data(:,8);      
elseif use_data ==  2        % IPI      
    ylev_data   =   csvread('Data_File_Output.csv');
    ylev        =   log(ylev_data(:,2));  
elseif use_data ==  3         % Hours
    ylev_data   =   csvread('Data_File_Output.csv');
    ylev        =   (ylev_data(:,6));  
elseif use_data ==  4          % Unemployment
    ylev_data   =   csvread('Data_File_Output.csv');
    ylev        =   log(ylev_data(:,8));  
end

% Defining the sample:
st              =   148;                    % = 1 for full sample 
                                            % = 148 for 1984Q1 
                                            % = 52 for 1960Q1. 
                
%% Options for local projections.

h1              =   1;      % For splines.        
lambda          =   100;    % Value of penalty of SLP.
r               =   3;      % r-1 = order of the limit polynomial for the SLP.      

%% Creating variables.

% LHS:
ylev2           =   ylev(st:end);
% Productivity:
zlev2           =   zlev(st:end);
% Inflation:
pilev2          =   pilev(st:end);
% Dummy for ZLB:
zlb2            =   zlb(st:end);
% Time trend and quadratic time trend: 
time2           =   time(st:end);
time222         =   time22(st:end);

% Creating W matrix:
ff              =   [ylev2 zlev2 pilev2];
w               =   lagmatrix(ff,1:P) ;     % Remove period t observations.

% Creating lagged variables interacted with period-t ZLB dummy:
[TY,YT]         =   size(w);
zlb3            =   zeros(TY,YT);
for jx = 1:YT
    zlb3(:,jx) = zlb2;
end
www             =   zlb3.*w; 
zlb4            =   1- zlb3;
wwww            =   zlb4.*w;

%% Creating the IRFs.

% LHS variable:   
y                       =   ylev2(1+P:end,:); 

% Local projection.
        % Standard case. 
w                       =   [zlb2(1+P:end) time2(1+P:end) time222(1+P:end)...
                             zlb2(1+P:end).*zlev2(1+P:end), wwww(1+P:end,:), ...
                             www(1+P:end,:)]; 
x                       =   (1-zlb2(1+P:end)).*zlev2(1+P:end,:); % Shock variable. 
y                       =   ylev2(1+P:end,:);  
lp                      =   locproj_var(y,x,w,h1,H,'reg'); 
% Extract the IRF from local projection: 
irf                     =   lp.IR(2:H+1); 
% Construct standard errors:
u                       =   lp.Y - lp.X*lp.theta;
[nwse_lp,vc_lp]         =   NeweyWest(u,lp.X,0,0);
        % ZLB case.
w                       =   [zlb2(1+P:end) (1-zlb2(1+P:end)).*zlev2(1+P:end,:)...
                             time2(1+P:end) time222(1+P:end), wwww(1+P:end,:),...
                             www(1+P:end,:)];
x                       =   zlb2(1+P:end).*zlev2(1+P:end);
lp_zlb                  =   locproj_var(y,x,w,h1,H,'reg'); 
% Extract the IRF from local projection: 
irf_zlb                 =   lp_zlb.IR(2:H+1); 
% Construct standard errors:
u                       =   lp_zlb.Y - lp_zlb.X*lp_zlb.theta;
[nwse_lp_zlb,vc_lp_zlb] =   NeweyWest(u,lp_zlb.X,0,0);

% Smooth local projection.
        % Standard case.  
w                       =   [zlb2(1+P:end) zlb2(1+P:end).*zlev2(1+P:end),...
                             time2(1+P:end) time222(1+P:end) wwww(1+P:end,:),...
                             www(1+P:end,:)]; 
x                       =   (1-zlb2(1+P:end)).*zlev2(1+P:end,:); % Shock variable. 
slp                     =   locproj_var(y,x,w,h1,H,'smooth',r,lambda);
% Extract the IRF from local projection: 
irf_slp                 =   slp.IR(2:H+1);
P2                      =   slp.P2;
% Construct standard errors:
u                       =   slp.Y - slp.X*slp.theta;
[nwse_slp,vc_slp,Q]     =   NeweyWest_slp(u,slp.X,0,0);
V2                      =   slp.B*vc_slp(1:slp.K,1:slp.K)*slp.B';
se_slp                  =   sqrt(diag(V2));
        % ZLB case.
w                       =   [zlb2(1+P:end) (1-zlb2(1+P:end)).*zlev2(1+P:end),...
                             time2(1+P:end) time222(1+P:end) wwww(1+P:end,:),...
                             www(1+P:end,:)];
x                       =   zlev2(1+P:end,:).*zlb2(1+P:end); % Shock variable. 
slp_zlb                 =   locproj_var(y,x,w,h1,H,'smooth',r,lambda); 
% Extract the IRF from local projection: 
irf_slp_zlb             =   slp_zlb.IR(2:H+1);
% Construct standard errors:
u                       =   slp_zlb.Y - slp_zlb.X*slp_zlb.theta;
[nwse_slp_zlb,vc_slp_zlb] = NeweyWest_slp(u,slp_zlb.X,0,0);
V22                     =   slp_zlb.B*vc_slp_zlb(1:slp_zlb.K,1:slp_zlb.K)*slp_zlb.B';
se_slp_zlb              =   sqrt(diag(V22));

%% Saving results to Excel.
p                       =   0.9;
bb                      =   tinv(p,TY-P);

irf_lp_baseline_LB      =   (irf-bb*nwse_lp(1:H));
irf_lp_baseline_UB      =   (irf+bb*nwse_lp(1:H));
irf_slp_baseline_LB     =   (irf_slp-bb*se_slp(1:H));
irf_slp_baseline_UB     =   (irf_slp+bb*se_slp(1:H));
irf_lp_zlb_baseline_LB  =   (irf_zlb-bb*nwse_lp_zlb(1:H));
irf_lp_zlb_baseline_UB  =   (irf_zlb+bb*nwse_lp_zlb(1:H));
irf_slp_zlb_baseline_LB =   (irf_slp_zlb-bb*se_slp_zlb(1:H));
irf_slp_zlb_baseline_UB =   (irf_slp_zlb+bb*se_slp_zlb(1:H));

xlswrite('Robustness.xlsx',irf,'Time Quadratic','B2');
xlswrite('Robustness.xlsx',irf_zlb,'Time Quadratic ZLB','B2');
xlswrite('Robustness.xlsx',irf_slp,'Time Quadratic','C2');
xlswrite('Robustness.xlsx',irf_slp_zlb,'Time Quadratic ZLB','C2');

%% Creating graph. 
     
t                       =   0:H-1;

if graph == 1
    figure
    plot(t,irf,'-k',t,irf_slp,'-b',t,irf_zlb,':k',t,irf_slp_zlb,':b','Linewidth',2)
    grid on 
    h = legend('LP','SLP','ZLB','ZLB SLP','Location','Best');
    set(h,'Interpreter','latex');
    if use_data == 1
        title('Output','Interpreter','latex')
    end
    if use_data == 2
        title('IP','Interpreter','latex')
    end   
    if use_data == 3
        title('Hours','Interpreter','latex')
    end 
    if use_data == 4
        title('Unemployment','Interpreter','latex')
    end 
    xlabel('$h$','Interpreter','latex')

    figure
    plot(t,(1/100)*irf_slp,'-b',t,(1/100)*irf_slp_zlb,':b','Linewidth',2)
    grid on 
    h = legend('SLP','ZLB SLP','Location','Best');
    set(h,'Interpreter','latex');
    if use_data == 1
        title('Output','Interpreter','latex')
    end
    if use_data == 2
        title('IP','Interpreter','latex')
    end   
    if use_data == 3
        title('Hours','Interpreter','latex')
    end 
    if use_data == 4
        title('Unemployment','Interpreter','latex')
    end 
    xlabel('$h$','Interpreter','latex')
end 